require(cregg)
require(RColorBrewer)
Set1 <- brewer.pal(9, "Set1")
windowsFonts(Helvetica = windowsFont("Helvetica"))

attribute.label <- c("Parents", "Identity", "Activity", 
                     "Policy", "Habit", "Intention")
level.label <- list(c("Independent", "Same party", "Opponent party"), 
                    c("Does not matter", "Personal insult"), 
                    c("Never", "Rarely", "Sometimes", "Often", "Very often"), 
                    c("Does not approve", "Partially approves", "Completely approves"), 
                    c("Always abstained", "Sometimes abstained", 
                      "Sometimes deviated", "Always loyal"), 
                    c("Abstain", "Same party", "Opponent party"))

choice.formula <- formula(choice ~ a.parents + b.identity + c.activity + d.policy + e.habit + f.intention)
rating.formula <- formula(rating ~ a.parents + b.identity + c.activity + d.policy + e.habit + f.intention)

#### load dataset ####
load("clean_data.Rdata")

# exclude Independents
partisan.data <- subset(clean.data, PID != 3)

#### analysis using marginal means (entire sample) ####
MM.choice.total.raw.result <- mm(partisan.data, choice.formula, 
                                 id = ~ respondent.id)
MM.choice.total.result <- MM.choice.total.raw.result
MM.choice.total.result$estimate <- MM.choice.total.result$estimate - mean(partisan.data$choice)
MM.choice.total.result$lower <- MM.choice.total.result$lower - mean(partisan.data$choice)
MM.choice.total.result$upper <- MM.choice.total.result$upper - mean(partisan.data$choice)

MM.rating.total.raw.result <- mm(partisan.data, rating.formula, 
                                 id = ~ respondent.id)
MM.rating.total.result <- MM.rating.total.raw.result
MM.rating.total.result$estimate <- MM.rating.total.result$estimate - mean(partisan.data$rating)
MM.rating.total.result$lower <- MM.rating.total.result$lower - mean(partisan.data$rating)
MM.rating.total.result$upper <- MM.rating.total.result$upper - mean(partisan.data$rating)

cairo_pdf("Figure_A2.pdf", width = 6.2, height = 3.8, pointsize = 7, family = "Helvetica")
layout(matrix(1:2, 1, 2))
par(mar = c(2, 0, 2, 0))
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 0.4), ylim = c(0, 26), 
     main = "Choice", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.4
loc.labels <- -1.07
loc.factors <- -1.2
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  for (j in 1:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(MM.choice.total.result$estimate[coefficients.location], plot.location, 
           pch = 19)
    segments(MM.choice.total.result$lower[coefficients.location], 
             plot.location, 
             MM.choice.total.result$upper[coefficients.location], 
             plot.location, lwd = 1)
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.4, -0.2, 0, 0.2, 0.4))
plot(NULL, NULL, bty = "n", xlim = c(-2.4, 0.8), ylim = c(0, 26), 
     main = "Rating", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.8
loc.labels <- -2.14
loc.factors <- -2.4
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  for (j in 1:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(MM.rating.total.result$estimate[coefficients.location], plot.location, 
           pch = 19)
    segments(MM.rating.total.result$lower[coefficients.location], 
             plot.location, 
             MM.rating.total.result$upper[coefficients.location], 
             plot.location, lwd = 1)
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.8, -0.4, 0, 0.4, 0.8))
dev.off()

#### analysis using marginal means (PID subgroups) ####
MM.choice.Dem.raw.result <- mm(subset(partisan.data, PID == 1), choice.formula, 
                               id = ~ respondent.id)
MM.choice.Dem.result <- MM.choice.Dem.raw.result
MM.choice.Dem.result$estimate <- MM.choice.Dem.result$estimate - mean(subset(partisan.data, PID == 1)$choice)
MM.choice.Dem.result$lower <- MM.choice.Dem.result$lower - mean(subset(partisan.data, PID == 1)$choice)
MM.choice.Dem.result$upper <- MM.choice.Dem.result$upper - mean(subset(partisan.data, PID == 1)$choice)

MM.rating.Dem.raw.result <- mm(subset(partisan.data, PID == 1), rating.formula, 
                               id = ~ respondent.id)
MM.rating.Dem.result <- MM.rating.Dem.raw.result
MM.rating.Dem.result$estimate <- MM.rating.Dem.result$estimate - mean(subset(partisan.data, PID == 1)$rating)
MM.rating.Dem.result$lower <- MM.rating.Dem.result$lower - mean(subset(partisan.data, PID == 1)$rating)
MM.rating.Dem.result$upper <- MM.rating.Dem.result$upper - mean(subset(partisan.data, PID == 1)$rating)

MM.choice.Rep.raw.result <- mm(subset(partisan.data, PID == 2), choice.formula, 
                               id = ~ respondent.id)
MM.choice.Rep.result <- MM.choice.Rep.raw.result
MM.choice.Rep.result$estimate <- MM.choice.Rep.result$estimate - mean(subset(partisan.data, PID == 2)$choice)
MM.choice.Rep.result$lower <- MM.choice.Rep.result$lower - mean(subset(partisan.data, PID == 2)$choice)
MM.choice.Rep.result$upper <- MM.choice.Rep.result$upper - mean(subset(partisan.data, PID == 2)$choice)

MM.rating.Rep.raw.result <- mm(subset(partisan.data, PID == 2), rating.formula, 
                               id = ~ respondent.id)
MM.rating.Rep.result <- MM.rating.Rep.raw.result
MM.rating.Rep.result$estimate <- MM.rating.Rep.result$estimate - mean(subset(partisan.data, PID == 2)$rating)
MM.rating.Rep.result$lower <- MM.rating.Rep.result$lower - mean(subset(partisan.data, PID == 2)$rating)
MM.rating.Rep.result$upper <- MM.rating.Rep.result$upper - mean(subset(partisan.data, PID == 2)$rating)

# Figure A.3
cairo_pdf("Figure_A3.pdf", width = 6.2, height = 3.8, pointsize = 7, family = "Helvetica")
layout(matrix(1:2, 1, 2))
par(mar = c(2, 0, 2, 0), xpd = TRUE)
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 0.4), ylim = c(0, 26), 
     main = "Choice", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.4
loc.labels <- -1.07
loc.factors <- -1.2
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  for (j in 1:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(MM.choice.Dem.result$estimate[coefficients.location], plot.location + 0.15, 
           pch = 19, col = Set1[2])
    segments(MM.choice.Dem.result$lower[coefficients.location], 
             plot.location + 0.15, 
             MM.choice.Dem.result$upper[coefficients.location], 
             plot.location + 0.15, col = Set1[2])
    points(MM.choice.Rep.result$estimate[coefficients.location], plot.location - 0.15, 
           pch = 15, col = Set1[1])
    segments(MM.choice.Rep.result$lower[coefficients.location], 
             plot.location - 0.15, 
             MM.choice.Rep.result$upper[coefficients.location], 
             plot.location - 0.15, col = Set1[1])
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.4, -0.2, 0, 0.2, 0.4))
legend(0.4, 26, legend = c("Self-identified D", "Self-identified R"), 
       pch = c(19, 15), col = Set1[c(2, 1)], bty = "n", xjust = 1, yjust = 0, cex = 0.8)
plot(NULL, NULL, bty = "n", xlim = c(-2.4, 0.8), ylim = c(0, 26), 
     main = "Rating", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.8
loc.labels <- -2.14
loc.factors <- -2.4
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  for (j in 1:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(MM.rating.Dem.result$estimate[coefficients.location], plot.location + 0.15, 
           pch = 19, col = Set1[2])
    segments(MM.rating.Dem.result$lower[coefficients.location], 
             plot.location + 0.15, 
             MM.rating.Dem.result$upper[coefficients.location], 
             plot.location + 0.15, col = Set1[2])
    points(MM.rating.Rep.result$estimate[coefficients.location], plot.location - 0.15, 
           pch = 15, col = Set1[1])
    segments(MM.rating.Rep.result$lower[coefficients.location], 
             plot.location - 0.15, 
             MM.rating.Rep.result$upper[coefficients.location], 
             plot.location - 0.15, col = Set1[1])
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.8, -0.4, 0, 0.4, 0.8))
legend(0.8, 26, legend = c("Self-identified D", "Self-identified R"), 
       pch = c(19, 15), col = Set1[c(2, 1)], bty = "n", xjust = 1, yjust = 0, cex = 0.8)
dev.off()